Se trabaja con datos reales obtenidos por la Intendencia de Montevideo, a través del Centro de Gestión de Movilidad, mediante sus contadores de autos ubicados en diferentes puntos de la ciudad. Estos datos vienen en intervalos acumulados de a 5 minutos y son captados por separado por cada carril que corre en un mismo sentido.
Los datos relevados haciendo uso de estos dispositivos, son de dominio público y la Intendencia los disponibiliza a través del sitio catálogo de datos que pertenece al estado.
El set de datos tiene muchas medidas con una frecuencia muy alta, muy útil y motivante su utilización.
Se pretende ajustar un modelo Poisson de media dependiente del tiempo para la cantidad de autos que transitan por una calle a lo largo de un día, obteniendo una distribución de probabilidad que nos permite simular estos eventos.
Originalmente los datos se presentan por archivos que acumulan períodos grandes de tiempo (meses, semestres, etc), variando según el año de relevamiento.
Cada uno de estos archivos posee, además del volumen de autos acumulados cada 5 minutos, el dato acumulado durante la última hora y otros parámetros que se muestran en la tabla:
| Dato | Tipo | Descripción | |-------|---------|-------------| | id_detector | Entero | Indica el id de la cámara que está monitoerando un determinado carril donde realiza el análisis de imagen para detectar la presencia de un vehículo | | id_carril | Entero | Es el número del carril que se está monitoreando. Va de 1 en adelante, siendo 1 el carril de más a la izquierda | | fecha | AAAA-MM-DD | Fecha en que fue tomada la muestra | | hora | hh:mm:ss | Hora en que fue tomada la muestra | | dsc_avenida | Texto | Nombre de la vía en la que se mide el tránsito | | dsc_int_anterior | Texto | Nombre de la vía que forma el cruce desde donde vienen los vehículos | | dsc_int_siguiente | Texto | Nombre de la vía que forma el cruce donde está el medidor, en general el mismo se encuentra un poco antes de esta vía. El sentido de circulación será desde el curce con dsc_int_anterior hacia el curce con dsc_int_siguiente. | | latitud | Float | latitud de donde está el medidor | | longitud | Float | longitud de donde está el medidor | | volumen | Enterno | Cantidad de vehículos detectados en el carril en los últimos 5 minutos | | volumen_hora | Enterno | Cantidad de vehículos detectados en el carril en la última hora |
Para el trabajo obligatorio se realizaron varias transformaciones previas sobre los datos:
Se tomó en cuenta la zona de Sarmiento y Av. José Requena y García con dirección a la rambla. Por lo que se tuvo que filtrar los archivos de varios meses, tomando solo los datos referidos a este punto.
Se agregan los datos de los 3 carriles allí existentes, contando únicamente los autos que pasaron en total.
Se tomó en cuenta la fecha y hora para lograr armar una serie de tiempo
Se completaron los datos faltantes con la media de la serie con las que se contaba originalmente.
# Instalaciones
#install.packages("data.table")
#install.packages("caret")
#Bibliotecas
library(astsa)
library(forecast)
library(data.table)
#Ajuste de tamaño de gráficas
options(repr.plot.width=15, repr.plot.height=8)
# Ejemplo de los datos crudos descargados de la web
df_ejemplo <- read.csv("./dataframe_ejemplo.csv", )
head(df_ejemplo)
| id_detector | id_carril | fecha | hora | dsc_avenida | dsc_int_anterior | dsc_int_siguiente | latitud | longitud | volumen | volumen_hora | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| <int> | <int> | <chr> | <chr> | <chr> | <chr> | <chr> | <dbl> | <dbl> | <int> | <int> | |
| 0 | 1 | 3 | 2017-01-01 | 00:00:00 | Sarmiento | Requena | Rambla | -34.91537 | -56.16836 | 0 | 0 |
| 1 | 1 | 3 | 2017-01-01 | 00:05:00 | Sarmiento | Requena | Rambla | -34.91537 | -56.16836 | 0 | 0 |
| 2 | 1 | 3 | 2017-01-01 | 00:10:00 | Sarmiento | Requena | Rambla | -34.91537 | -56.16836 | 6 | 72 |
| 3 | 1 | 3 | 2017-01-01 | 00:15:00 | Sarmiento | Requena | Rambla | -34.91537 | -56.16836 | 2 | 24 |
| 4 | 1 | 3 | 2017-01-01 | 00:20:00 | Sarmiento | Requena | Rambla | -34.91537 | -56.16836 | 3 | 36 |
| 5 | 1 | 3 | 2017-01-01 | 00:25:00 | Sarmiento | Requena | Rambla | -34.91537 | -56.16836 | 6 | 72 |
Se elige la avenida Sarmiento para esta tarea, donde el mismo cuenta con 3 detectores en total, cubriendo cada uno de sus carriles. Luego de un procesamiento, se llega a un dataframe con una columna de fechas date y otra con volumen. Este procesamiento se hace fuera de la Notebook con el script main.py, donde se saca la infomación de cada csv, en este caso, todo 2017.
# Datos crudos obtenidos cada 5 minutos de cada sensor de carril
df_raw = read.csv("./2017_2018_raw.csv", sep =";")
head(df_raw)
tail(df_raw)
| date | volumen | |
|---|---|---|
| <chr> | <int> | |
| 1 | 1/1/17 00:00 | 0 |
| 2 | 1/1/17 00:00 | 0 |
| 3 | 1/1/17 00:00 | 0 |
| 4 | 1/1/17 00:05 | 0 |
| 5 | 1/1/17 00:05 | 0 |
| 6 | 1/1/17 00:05 | 0 |
| date | volumen | |
|---|---|---|
| <chr> | <int> | |
| 307837 | 31/12/17 23:50 | 0 |
| 307838 | 31/12/17 23:50 | 3 |
| 307839 | 31/12/17 23:50 | 0 |
| 307840 | 31/12/17 23:55 | 0 |
| 307841 | 31/12/17 23:55 | 3 |
| 307842 | 31/12/17 23:55 | 0 |
freq_5min = 365 * 24 * 60 / 5
freq_5min
freq_15min = 365 * 24 * 60 / 15
freq_15min
freq_sem_15min = 365 * 24 * 60 / ( 3 * 5 * (365 / 7) )
freq_sem_15min
Se convierte a una serie de tiempo los datos cargados, sumarizándolos cada 15 minutos.
# Generar una serie de tiempo con datos cada 5 minutos
df_ts_raw = ts(df_raw$volumen, frequency=freq_5min, start=c(2017,1))
# Agregar la serie anterior cada 15 minutos (sumando los volúmenes)
df_aggregate_raw = aggregate(df_ts_raw, nfrequency = freq_15min, FUN = sum)
df_aggregate_ts_raw = ts(df_aggregate_raw, frequency = freq_15min*3, start=c(2017,1))
# Plot de la TS
tsplot(df_aggregate_ts_raw, main="Serie original sumarizada cada 15 minutos")
En esta serie se puede ver que no alcanza el 2018, porque faltan datos en varios momentos del año, para corregir esto, se generó una serie anual con la frecuencia deseada y se hico un merge con los datos originales.
# Copia de dataframe original
originaldf <- df_raw
# Creación de nueva columna para new date
originaldf$minAsPOSIX <- as.POSIXct(originaldf$date, format="%d/%m/%y %H:%M",tz="GMT")
# Generación del vector para cada 5 minutos
ndays <- 365 # días a generar
minAsNumeric <-seq(0, 60 * 60 * 24 * ndays, by = 60 * 5) # de 0 a 365 días con un paso de 5 minutos
# Conversión de secuencia a posix (formato calendario)
minAsPOSIX <- as.POSIXct(minAsNumeric, origin = "2017-01-01", tz = "GMT")
# Columna con datos calendario para 1 año
base_df = data.frame(minAsPOSIX)
# Merge de datos calendario completos con los datos de volumen original
df_merged <- merge(base_df, originaldf, all.x = TRUE, by="minAsPOSIX")
# Agregación de datos de volumen por los carriles (los que tengan misma fecha y hora los suma)
newdf = aggregate(df_merged$volumen, by = list(Date=df_merged$minAsPOSIX), FUN=sum)
colnames(newdf) <- c('date', 'volumen')
ts = ts(newdf$volumen, frequency=freq_5min, start=c(2017,1))
tsplot(ts, main="Serie original con NAs donde faltan datos")
Ahora se obtiene un nuevo dataframe donde en la primer columna están los datos *date* de fecha completos en todo el tiempo deseado sin huecos, con su correspondiente sumarización en la columna *volumen*.
En caso de faltar datos, quedan NA, luego serán sustituidos por la media de la serie a trabajar.
Se analiza en períodos más cortos para visualizar la forma, reconocer mejor los días, las semanas y los fines de semana.
NOTA: El 1/1/2017 fue domingo.
# Muestras por día (1 muestra cada 5 minutos x 2)
samples_per_day = (60 / 5) * 24
samples_per_week = samples_per_day * 7
# Días a tomar para análisis más pequeño
days_to_crop = 21
initial_week = 35
# Intervalo
ini = initial_week * samples_per_week
fin = ini + samples_per_day * days_to_crop
st = 2017 + ini/length(ts)
# Serie recortada
df_crop = ts[ini : fin]
# Promedio de serie cortada y sustituir
mean = mean(df_crop, na.rm=TRUE)
print(paste("Media:", mean))
df_crop[is.na(df_crop)] = mean
df_ts_crop = ts(df_crop, start = st, freq = freq_5min)
[1] "Media: 29.6276912885061"
# Plot de la serie cortada
tsplot(df_ts_crop, main="Serie con marcas semanales")
# Plot barras en las semanas
lp = 2
picos = rep(c( 0, rep(150, lp/2), 0 ,rep(NA, samples_per_week - lp - 4), 0 , rep(150, lp/2), 0 ), days_to_crop/7)
ts_picos = ts(picos, start = st, freq = freq_5min)
lines(ts_picos, lwd = 3, col = 4)
Visto que la serie difiere mucho con las estaciones del año, se decide utilizar el último 30% de la serie (setiembre a diciembre aproximadamente).
# Intervalo
ini = initial_week * samples_per_week
fin = length(ts)
st = 2017 + ini/length(ts)
# Serie recortada
df_crop = ts[ini : fin]
# Media de serie cortada y sustituir los NAs
mean = mean(df_crop, na.rm=TRUE)
print(paste("Media:", mean))
df_crop[is.na(df_crop)] = mean
ts_crop = ts(df_crop, start = st, freq = freq_5min)
[1] "Media: 30.2471385934919"
# Plot de la serie cortada
tsplot(ts_crop, main="Serie original cortada con marcas semanales")
# Plot barras en las semanas
lp = 2
picos = rep(c( 0, rep(150, lp/2), 0 ,rep(NA, samples_per_week - lp - 4), 0 , rep(150, lp/2), 0 ), 52 - initial_week)
ts_picos = ts(picos, start = st, freq = freq_5min)
lines(ts_picos, lwd = 3, col = 4)
acf(ts_crop)
pacf(ts_crop)
Del ACF se puede observar que hay una correlación alta con valores pasados, pero el PACF nos dice cuánto para atrás deberíamos de ir para ajustar la serie.
# Estudio de la periodicidad
samples_per_day = 288 # 24 * 60 / 5
s = ts_crop
s = s-mean(s) # Se quita la media
n = length(s)
I = abs(fft(s))^2 # FFT y módulo cuadrado
I = I[1:floor(n/2)] # Recorto el vector a las frecuencias observables
P = (4/n^2)*I # Escalado del periodograma
f = (0:(n/2-1))/n*(freq_5min) # Vector de frecuencias para hacer el gráfico (se multiplica por la frecuencia de la serie, para normalizar
# Ploteo del periodograma
plot(f, P, type="b", xlab="Frecuencia", main="Periodograma escalado", col=4, lwd=1, pch=19)
plot(f[0:1000], P[0:1000], type="b", xlab="Frecuencia", main="Periodograma escalado", col=4, lwd=1, pch=19)
Del periodograma se pueden extraer las frecuencias principales. Buscamos los máximos del vector P para obtener el valor de las frecuencias.
samples_per_day = 1
# Primera frecuencia de mayor importancia
f1 = f[which.max(P)]*samples_per_day
print(paste("f1:",f1))
print(paste("which.max(P1):",which.max(P)))
# Segunda frecuencia de mayor importancia
P[which.max(P)] = 0
f2 = f[which.max(P)]*samples_per_day
print(paste("f2:",f2))
print(paste("which.max(P2):",which.max(P)))
# Tercera frecuencia de mayor importancia
P[which.max(P)] = 0
f3 = f[which.max(P)]*samples_per_day
print(paste("f3:",f3))
print(paste("which.max(P3):",which.max(P)))
# Cuarta frecuencia de mayor importancia
P[which.max(P)] = 0
f4 = f[which.max(P)]*samples_per_day
print(paste("f4:",f4))
print(paste("which.max(P4):",which.max(P)))
# Quinta frecuencia de mayor importancia
P[which.max(P)] = 0
f5 = f[which.max(P)]*samples_per_day
print(paste("f5:",f5))
print(paste("which.max(P5):",which.max(P)))
[1] "f1: 364.978878537122" [1] "which.max(P1): 121" [1] "f2: 729.957757074243" [1] "which.max(P2): 241" [1] "f3: 1094.93663561137" [1] "which.max(P3): 361" [1] "f4: 416.684219663214" [1] "which.max(P4): 138" [1] "f5: 313.273537411029" [1] "which.max(P5): 104"
t = time(ts_crop)
t_start_crop = t[1]
#f = 365
# Creación de fit con las 2 frecuencias mas significativas halladas
fit_df_ts_crop = lm(ts_crop ~ cos(2*pi*f1*t)+sin(2*pi*f1*t)+cos(2*pi*f2*t)+sin(2*pi*f2*t)+cos(2*pi*f3*t)+sin(2*pi*f3*t))
#fit_df_ts_crop = lm(ts_crop ~ cos(2*pi*f*t)+sin(2*pi*f*t)+cos(4*pi*f*t)+sin(4*pi*f*t)+cos(6*pi*f*t)+sin(6*pi*f*t)) # en caso de mantener f fija
summary(fit_df_ts_crop)
Call:
lm(formula = ts_crop ~ cos(2 * pi * f1 * t) + sin(2 * pi * f1 *
t) + cos(2 * pi * f2 * t) + sin(2 * pi * f2 * t) + cos(2 *
pi * f3 * t) + sin(2 * pi * f3 * t))
Residuals:
Min 1Q Median 3Q Max
-62.319 -7.241 -0.529 7.595 76.616
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 30.24714 0.08896 340.02 <2e-16 ***
cos(2 * pi * f1 * t) 20.01699 0.12580 159.11 <2e-16 ***
sin(2 * pi * f1 * t) -8.80723 0.12580 -70.01 <2e-16 ***
cos(2 * pi * f2 * t) -10.81332 0.12580 -85.95 <2e-16 ***
sin(2 * pi * f2 * t) -2.34551 0.12580 -18.64 <2e-16 ***
cos(2 * pi * f3 * t) -3.85440 0.12580 -30.64 <2e-16 ***
sin(2 * pi * f3 * t) 7.18327 0.12580 57.10 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 16.54 on 34555 degrees of freedom
Multiple R-squared: 0.5495, Adjusted R-squared: 0.5494
F-statistic: 7025 on 6 and 34555 DF, p-value: < 2.2e-16
De el ajuste se puede saber que todas las componentes son significativas.
adjustment_fit_df_ts_crop = ts(fitted(fit_df_ts_crop), start=t_start_crop, freq=freq_5min)
tsplot(ts_crop, col=4, lwd=2, main = "Serie original vs ajuste lineal")
lines(adjustment_fit_df_ts_crop, lwd=2)
# Residuals del ajuste
residuals_fit = residuals(fit_df_ts_crop)
samples_per_day = 288
# Ajuste para las primeras 2 semanas del dataframe recortado
tsplot(head(ts_crop, samples_per_day*days_to_crop), col=4, lwd=2, main = "Principio de serie original vs ajuste lineal aumentado")
lines(head(adjustment_fit_df_ts_crop, samples_per_day*days_to_crop), lwd=2)
# Ajuste para las últimas 2 semanas del dataframe
tsplot(tail(ts_crop, samples_per_day*days_to_crop), col=4, lwd=2, main = "Final de serie original vs ajuste lineal aumentado")
lines(tail(adjustment_fit_df_ts_crop, samples_per_day*days_to_crop), lwd=2)
acf(residuals_fit)
pacf(residuals_fit)
checkresiduals(residuals_fit)
Warning message in modeldf.default(object): “Could not find appropriate degrees of freedom for this model.”
Sigue habiendo una componente fuerte y 2 menores a pesar de haber ajustado la serie.
# Auto ARIMA con Seasonal
autoarima = auto.arima(ts_crop, seasonal = TRUE)
ajuste_autoarima = ts(fitted(autoarima), start=t_start_crop, freq=freq_5min)
summary(autoarima)
residuals_autoarima = residuals(autoarima)
Series: ts_crop
ARIMA(5,0,0) with non-zero mean
Coefficients:
ar1 ar2 ar3 ar4 ar5 mean
0.4533 0.2501 0.1449 0.0760 0.0514 30.2471
s.e. 0.0054 0.0059 0.0060 0.0059 0.0054 1.6497
sigma^2 = 56.14: log likelihood = -118643.6
AIC=237301.3 AICc=237301.3 BIC=237360.5
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.0002061926 7.491732 5.370428 -Inf Inf NaN -0.0001746736
acf(residuals_autoarima, na.action = na.pass)
pacf(residuals_autoarima, na.action = na.pass)
Con el Auto ARIMA se ve como el ACF y el PACF quedan limpios, sin encontrar correlación con valores pasados luego de aplicar un arima autoregresivo de orden 5, media móvil 0 e integrador 0.
Se quiere armar una lista con vectores que representan el conteo de autos para determinado momento.
# Lista de vectores
new_daily_vector = list()
# Muestras por día
samples_per_day_5min = 288
# Frecuencia a contemplar para el muestreo
days_to_sum = 7
# Muestras totales por días usados
samples_per_sum = samples_per_day_5min * days_to_sum
# Cuántos días se quiere muestrear
data_forward = 12
for (i in 1:data_forward){
if (i != 1){
new_daily_vector[[i]] = ts_crop[(i-1)*(samples_per_sum+1)]
}else{
new_daily_vector[[i]] = ts_crop[1]
}
for (j in 1:samples_per_sum){
new_daily_vector[[i]] = append(new_daily_vector[[i]], round(ts_crop[(i-1)*samples_per_sum + j]),0) # Redondeo sino optim da Warning
}
new_daily_vector[[i]] = rev(new_daily_vector[[i]])
}
# Plot vectores diarios en misma grafica
ts1 = ts(new_daily_vector[[1]], start=0, freq=samples_per_day_5min)
tsplot(ts1, ylim = c(-1,130), main = "Superposición de días a modelar")
# Plot de los vectores hallados
for (i in 2:data_forward){
ts = ts(new_daily_vector[[i]], start=0, freq=samples_per_day_5min)
lines(ts, col = i)
}
# Transponer la matríz para tener en las filas las muestras para el mismo momento para la frecuencia determinada.
transposed_list = do.call(Map,c(f=c,new_daily_vector))
La distribución de Poisson se usa para modelar el número de eventos que ocurren en un proceso de Poisson. Sea $X∼P(λ)$, esto es, una variable aleatoria con distribución de Poisson donde el número medio de eventos que ocurren en un determinado intervalo con esta función de probabilidad:
$$P(X=x)= \frac{e^{-\lambda} \lambda^x}{x!}$$Esto se obtiene con la función *dpois*. A esta se le ingresa la serie a ajustar, la función lambda a utilizar, y si se devuelve la probabilidad como su logaritmo. Al ser eventos independientes, las probabilidades se multiplican. Como se obtiene el logaritmo, éstas se terminan sumando y haciéndose negativas para hallar el mínimo ya que se quiere maximizar la verosimilitud.
# Lista para guardar los parámetros de la salida optim
est_vector_par = list()
# Lista para guardar los valores medios iterados
mean_vector_par = list()
for (i in 1:samples_per_sum){
f <- function(lambda) {
-sum(dpois(unlist(transposed_list[i]),lambda = lambda, log = TRUE), na.rm = TRUE)
}
est = optim(par = transposed_list[i], fn = f, method = "Brent", lower = 0, upper = 150)
est_vector_par = append(est_vector_par, est$par)
mean_vector_par = append(mean_vector_par, mean(as.numeric(unlist(transposed_list[i]))))
}
est_vector_par_ts = ts(est_vector_par, frequency = samples_per_day_5min, start=c(0,1))
mean_vector_par_ts = ts(mean_vector_par, frequency = samples_per_day_5min, start=c(0,1))
# Plot de vectores hallados
tsplot(est_vector_par_ts, col = 4, main = "Resultado de optimizador para cada punto")
lines(mean_vector_par_ts, col = 2)
El resultado de la optimización es la media de la cantidad de autos en un mismo tiempo para cada ventana semanal. Lo cual tiene sentido, porque el lambda es eso, la media de ocurrencias del evento.
# Ploteo de Lambda encontrado frente a los datos analizados
tsplot(est_vector_par_ts, col = 4, lwd = 5, ylim = c(-1,130), main = "Superposición vs evaluación para cada punto")
for (i in 2:data_forward){
ts = ts(new_daily_vector[[i]], start=0, freq=samples_per_day_5min)
lines(ts)
}
lines(est_vector_par_ts, col = 4, lwd = 2)
Ahora se quiere encontrar una función sinusoidal que pueda estimar la serie.
Se busca estimar, a través de la optimización simulando una distribución de Poisson, la función $\lambda(t)$ que "mejor" ajuste la media de cada evento para obtener dicha serie.
Se quieren encontrar los parámetros:
$$Params = \{a, b, c, d, e\}$$Para la función:
$$\lambda(t) = a + b\ cos(2\pi f t) + c\ sin(2\pi f t) + d\ cos(4\pi f t) + e\ sin(4\pi f t)$$Con $f = 365$ en este caso para no agregar complejidad a la función.
Definimos el vector de parámetros iniciales aproximados a los antes encontrados por el modelado lineal *lm* para darle menos trabajo al optimizador.
init.par = c(30, -20, -8, 4, -7, 3, 7)#,0)
x = round(ts_crop)
t = time(ts_crop)
f <- function(params) {
lambdat = params[1] + params[2]*cos(2*pi*365*t) + params[3]*sin(2*pi*365*t) + params[4]*cos(4*pi*365*t) + params[5]*sin(4*pi*365*t)+ params[6]*cos(6*pi*365*t) + params[7]*sin(6*pi*365*t)
#lambdat = params[1] + params[2]*cos(2*pi*365*(t+params[8])) + params[3]*sin(2*pi*365*(t+params[8])) + params[4]*cos(4*pi*365*(t+params[8])) + params[5]*sin(4*pi*365*(t+params[8]))+ params[6]*cos(6*pi*365*(t+params[8])) + params[7]*sin(6*pi*365*(t+params[8]))
lambdat[which(lambdat < 0)] = 0
-sum(dpois(x, lambda = lambdat, log = TRUE))
}
est = optim(par = init.par, fn = f, method="BFGS", hessian=TRUE, control=list(trace=1, maxit=100))
est
initial value 229029.552541 iter 10 value 221087.821826 final value 221083.133073 converged
| 2345.29123 | 917.55915 | 952.31453 | -296.30262 | 914.3812 | -726.83419 | -53.95972 |
| 917.55915 | 1024.49427 | 457.19058 | 95.36251 | 449.1774 | -319.25510 | 139.54150 |
| 952.31453 | 457.19058 | 1320.79690 | -503.13711 | 822.1967 | -774.83971 | 22.95248 |
| -296.30262 | 95.36251 | -503.13711 | 1001.54179 | -317.6491 | 644.85112 | 207.24128 |
| 914.38121 | 449.17738 | 822.19666 | -317.64909 | 1343.7494 | -745.07324 | 272.70803 |
| -726.83419 | -319.25510 | -774.83971 | 644.85112 | -745.0732 | 1413.67222 | 26.17709 |
| -53.95972 | 139.54150 | 22.95248 | 207.24128 | 272.7080 | 26.17709 | 931.61898 |
Luego de algunas iteraciones, se encuentran los valores de los parámetros de ajuste. A simple viste se ve una diferencia leve en amplitudes halladas respecto a la LM. Se grafica la serie $Lambda$ encontrada contra la serie original.
ajuste_f <- function(t) {
lambdat = est$par[1] + est$par[2]*cos(2*pi*365*t) + est$par[3]*sin(2*pi*365*t) + est$par[4]*cos(4*pi*365*t) + est$par[5]*sin(4*pi*365*t) + + est$par[6]*cos(6*pi*365*t) + est$par[7]*sin(6*pi*365*t)
#lambdat = est$par[1] + est$par[2]*cos(2*pi*365*(t+est$par[8])) + est$par[3]*sin(2*pi*365*(t+est$par[8])) + est$par[4]*cos(4*pi*365*(t+est$par[8])) + est$par[5]*sin(4*pi*365*(t+est$par[8])) + + est$par[6]*cos(6*pi*365*(t+est$par[8])) + est$par[7]*sin(6*pi*365*(t+est$par[8]))
}
eval_poisson = eval(ajuste_f(t))
tsplot(ts_crop, col=4, lwd=2, main = "Serie original truncada vs evaluación de Lambda(t)")
lines(eval_poisson, col = 2)
tsplot(head(ts_crop,288*7), col=4, lwd=2, main = "Zoom de serie original truncada vs evaluación de Lambda(t)")
lines(eval_poisson, col = 2, lwd = 3)
# Residuales encontrados por la resta de la serie original contra la evacluacion de Lambda
residuals_poisson = ts_crop - eval_poisson
tsplot(ts(residuals_fit, start = 0, freq=freq_5min), col = 4, main = "Residuales de la LM vs Residuales en Evaluación de Lambda(t)")
lines(ts(residuals_poisson, start = 0, freq=freq_5min), col = 2)
tsplot(head(ts(residuals_fit, start = 0, freq=freq_5min), 288*7), col = 4, main = "Head de Residuales de la LM vs Residuales en Evaluación de Lambda(t)")
lines(head(ts(residuals_poisson, start = 0, freq=freq_5min), 288*7), col = 2)
acf(residuals_poisson)
pacf(residuals_poisson)
# Varianza para cada Residuals
print("Varianza de Residuals Poisson")
var(residuals_poisson)
print("Varianza de ajuste lineal")
var(residuals_fit)
[1] "Varianza de Residuals Poisson"
[1] "Varianza de ajuste lineal"
# Root Mean Square Error para cada Residuals
print("RMSE de Residuals Poisson")
sqrt(mean(residuals_poisson^2))
print("RMSE de ajuste lineal")
sqrt(mean(residuals_fit^2))
[1] "RMSE de Residuals Poisson"
[1] "RMSE de ajuste lineal"
Las Varianzas y RMSE dan similares, un poco menores en el LM.
Cuando se incluyó una fase dentro de las componentes senoidales, se obtuvo un parámetro en la estimación muy chico, que a la hora de evaluar la serie la Varianza dió 282, un poco peor que sin este desfazaje. Esta prueba se hizo para analizar si, al agregar esta componente no lineal, se lograban mejores resutados.
Se grafican a continuación ambos modelados para compararlos visualmente.
tsplot(ts(fitted(fit_df_ts_crop), start = t_start_crop, freq = freq_5min), main = "Gráfico de la LM y la Evacluación de Lambda")
lines(eval_poisson, col = 2)
one_week_ts_crop = head(ts_crop,samples_per_day_5min*days_to_sum)
tsplot(one_week_ts_crop, main = "Gráfico de la LM y la Evacluación de Lambda")
lines(head(ts(fitted(fit_df_ts_crop), start = t_start_crop, freq = freq_5min),288*7), col = 4, lwd = 4)
lines(head(eval_poisson,samples_per_day_5min*days_to_sum), col = 2, lwd = 4)
legend(x = min(time(one_week_ts_crop), time(one_week_ts_crop)), y = max(one_week_ts_crop,one_week_ts_crop), legend=c("Serie original", "Linear Model", "Evaluación Poisson"), fill = c(1,4,2))
Ambos modelados son muy similares, dado que los parámetros hallados también los son.
Se ajustan los residuales de la función Lambda a un autoarima.
# Auto ARIMA para los residuals de la evalcuación de Poisson
autoarima_poisson = auto.arima(residuals_poisson)
autoarima_poisson_ts = ts(fitted(autoarima_poisson), start=t_start_crop, freq=freq_5min)
summary(autoarima_poisson)
residuals_autoarima = residuals(autoarima_poisson)
Series: residuals_poisson
ARIMA(5,0,2) with non-zero mean
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2 mean
1.3683 -0.8478 0.2187 0.0788 0.1403 -0.9452 0.6921 -0.0007
s.e. 0.0260 0.0485 0.0185 0.0093 0.0115 0.0249 0.0408 0.7064
sigma^2 = 54.14: log likelihood = -118018.3
AIC=236054.6 AICc=236054.6 BIC=236130.7
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 5.849429e-05 7.357473 5.292776 -12.5111 243.4675 NaN -0.00541282
acf(residuals_autoarima)
pacf(residuals_autoarima)
checkresiduals(residuals_autoarima)
Warning message in modeldf.default(object): “Could not find appropriate degrees of freedom for this model.”
